THE QUADRATIC CHARACTER EXPERIMENT 

JEFFREY STOPPLE 

Abstract. A fast new algorithm is used compute the zeros of 
the quadratic character L-functions for all negative fundamental 
discriminants with absolute value 10^^ < d < 10^^ + 10''. These 
are compared to the 1-level density, including various lower order 
terms. These terms come from, on the one hand the Explicit For- 
mula, and on the other the L-functions Ratios Conjecture. The 
latter give a much better fit to the data, providing numerical evi- 
dence for the conjecture. 



1. Introduction. 

Predictions. Standard conjectures [5j predict that the low lying zeros 
of quadratic Dirichlet L-functions should be distributed according to a 
symplectic random matrix model. To make this more precise, we'll in- 
troduce some notation. Let Xd be a real, primitive character modulo d, 
and suppose furthermore that Xd{—^) -d is a fundamental discriminant. 
Let g{T) be a Schwartz class test function. Then the 1-level density for 
the zeros 1/2 -|- of L(s, Xd) should satisfy 

d<X 7d ^ ^ 

where X* is the cardinality of fundamental discriminants Xd{—^)'d with 
d < X. This is a theorem [8] if the support of the Fourier transform of 
g is suitably restricted. 

Recently Conrey and Snaith [2] made a precise prediction for the 
lower order arithmetic terms in the 1-level density. Their prediction is 
conditional, assuming the L-functions Ratios Conjecture [1]. Miller [7] 
then proved (under typical restrictions for the test function g^r)) that 
these lower order terms exist and agree with the prediction in [2]. 
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Experiments. Zeros of Dirichlet L-functions were first computed by 
Davies and Haselgrove, by Spira, and by Rumley. Rubinstein, as one 
portion of his thesis [10], was the first to compute enough low lying 
zeros to meaningfully test ([T]). However, the numerical methods devel- 
oped in [To] were optimized to compute L{l/2 + zt, Xd) for large real t 
rather than large d. 

This paper. The next section develops an algorithm to compute low 
lying zeros (0 < t < 1.) which is fast for large d. This is a modification 
of the idea behind [llj. The subsequent section has a discussion of 
the data from the computation of the zeros of approximately 3 • 10^ 
quadratic character L-functions for negative fundamental discriminants 
—d with d > 10.^^. This is followed by some implementation notes, and 
an Appendix on Miller's 'refined' 1-level density and the L-functions 
Ratios Conjecture. 

Acknowledgements. Thanks to Mike Rubinstein for sharing his data 
from pTOj, and David Farmer for pointing me towards [7] Thanks to 
both Steven J. Miller and to the anonymous referee for their careful 
reading of the manuscript and numerous helpful suggestions. 

2. Algorithm 

We are going to compute the L-function on the critical line by means 
of an approximate functional equation, an idea that goes back to Lavrik 
and was first implemented by Weinberger p!5] . With x a real character 
modulo d, let a = (1 — x(— 1))/2, and with t > use s to denote 
(1/2 + zt + a)/2Q Define 

(2) Z{t, x) = e(l/2 + zt, x) = Yl ^(^)^" 2 Re{G{s, nn^/d)), 

n 

where ^ 

x) = a;~*r(s, x) = ( exp(—yx)y''^ — . 

Ji y 

As in ^12], the tail of the series, the sum of terms n > A^, is bounded 
by d"^ exp{—N^TT /d) / (ttN)^ , so if we want to compute to D digits of 
accuracy, we should have 

(3) d^exp{-N\/d)/{7TN)'' KlO-"". 
Certainly 

(4) Ar>dl/2log(rf2^QD)l/2^-l/2 

^We are not actually assuming the Generalized Riemann Hypothesis, but we are 
only looking for zeros on the critical line. 
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would suffice; later we'll see we can do better given any particular d. 

Differentiating with respect to x under the integral defining G{s, x) 
we see that 

(5) ^G{s,x) = - j\M-xy)y'^'j = + l,x), 

while integration by parts, on the other hand, gives 

(6) G{s + 1, x) = !^Piz^ + '-G{s, x). 

X X 

Equations ([5]) and ([6]) give a nice recursive relation for all the derivatives 
x) in terms of This, in turn, motivates a consideration 

of Taylor expansions. 

Suppose we compute G{s, x) by a Taylor series expansion (in the 
second variable, centered at xq) to B terms, where 5 is a parameter to 
be determined. 

Lemma. We can bound the remainder in the Taylor expansion by a 
function Rb{x,xo) (defined below) which satisfies 



(7) Rb{x, Xo) < ^ ^, ' T{B, xo) < ^ 

Proof. We have 

/oo 
exp{~xy)y^dy, 

since s is in the critical strip. By the integral formula for the remainder 
in Taylor's theorem, we can bound that remainder by 
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Rb{x,xq) = — J J exp{-uy)y'^dy {x - u)^du. 
Change the order of integration and let t = x — m to get 

exp(— xy) / exp{—ty)t^dty^dy. 

J X — Xq 

Now integrate by parts in the t integral to get 

exp(-xy)(x - xo)'^exp((x - xo)y)y^~^dy 



-1 



1 

m 



- i?B_i(x,Xo) 
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Or, in other words, 

Rb{x, Xo) + Rb-i{x, Xq) = — G{B, xq) 

_ l) 

This imphes the first inequahty. For the second, we observe 



T{B,xo) = / exp{-y)y 

'xo y 



< rexpi-y)y^^=TiB) = iB-l)\ 

Jo y 



□ 

The first inequahty is stronger, so it's good for the actual computa- 
tion. The second is weaker, but simple enough to be useful in proving 
the theorem. 

Now we're ready to put the Taylor expansions to good use. Similar to 
the method of [11] , we partition the set {n^ 1 1 < < A^} into intervals 

for j = 1,...,T, where Fj is the jth Fibonacci number. We then 
compute the function G by a Taylor expansion in the second variable, 
centered at irFj/d, and truncated to B terms: 

B 

2Re(G(s,7™Vrf)) ^ ^G'j,fc(t)(7r/d)^ ■ {n^ - F.j)'' , 

k=0 

where 

(8) Gj-fc(t) = 2Re{G'^''\s,7TFj/d))/k\. 
Theorem 1. We can compute Z{t,x) as 

T B 

(9) Z(t,x) = 5^5^G„.(t)(vr/d)'= X{n)n\n'-F,f 

3=1 fe=0 n^e/j 

to D digits of accuracy, where T and B are both 0(log((i)), the implied 
constants depending on D. 

The expression 

(10) "=^- xinKin' - F,f 
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is a precomputation independent of s in integers which is 0{N ■ B) = 
0((i^/^ log((i)^). Subsequently, individual evaluations of Z{t,x) cost 
only 0(T- 5) = 0{log{d)^). 

Proof. Of course, the outermost sum on j <T and the innermost sum 
on n"^ in Ij combine to give the squares of all n < A^; the middle sum 
giving the needed Taylor expansions. We need A^^ to be in the last 
interval It, so 

< ^ $^+Vv^, where $ = i±^; 



with < d^/^+e ^]^jg implies that T < log(ci) suffices. 

This is all well and good, but we need to show that using Taylor ex- 
pansions at points spaced in what is essentially a geometric progression, 
does not require an unreasonable number of terms B in each expansion 
in order to compute accurately. Use < 1, < n, and the rough 

estimate cf^'^ for the L-series truncation parameter N. Assuming the 
errors we make in computing each G{s,'kv? / d) are independent with 
standard deviation e, then the standard error in the sum ^ is bounded 
by [3] 

1/2 

« e ■ d''\ 



2 



n=l 



where we approximated a sum by an integral. We want e ■ d^/'^ < 10 ^, 
or 

6 < lo-^^/-3/^ 

which will determine how many terms B we need in each Taylor ex- 
pansion. We'll use the weaker inequality in the Lemma with 

xq = TiFj/d, X < TcFj^i/d, 

which makes the error 

B 

- V 



Thus we want 
($ - 1)^ < 10-^ci-='/^ or 10^^=^/^ < ($ - 1)-^ 



B 



and so -B = 0{\og{d)) suffices. □ 
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Figure 1. Histogram of the lowest zero. 

For a single function evaluation (for example, determining whether 
Z{0.,x) > 0.) this algorithm is no improvement over [15]; summing 
the series requires 0{dlog{d)Y^'^ terms by (|4]). If one wants to do an 
arbitrarily large number of function evaluations, the improvement is 
spectacular: from exponential down to polynomial (in terms of the 
number of digits of d which is ~ \og{d)). This is deceptive, though, 
because what one really wants to do is find the all zeros with, say, 
< t < 1. (Larger t intervals requires computing G{s, x) via the meth- 
ods of [10] which in turn necessitates re-doing the precomputation.) 
Since there are 0{\og{d)) such zeros and each can be found with 0(1) 
evaluations, the precomputation still dominates as a theoretical result. 
But as Jan L. A. van de Snepscheutj^ wrote 

"In theory, there is no difference between theory and 
practice. But, in practice, there is." 

3. Data 

Zeros with t < 1. were computed for all of the 1,039,654 negative 
fundamental discriminants —d in the range 10^^ ^ d < 10^^ + 10^, a 
total of 12,202,567 zeros. Figure [l] shows a histogram for the imaginary 



^not Yogi Berra. 
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part of the lowest lying zero, rescaled by log(10^^)/(27r). The lowest 
zero found was at t = 0.0013104755 corresponding to the discriminant 
-1,000,008,582,815. 

Figure |2] shows the histogram of imaginary parts of all the zeros, 
again rescaled by log(10^^)/(27r) ~ 4.39761. The upper curve (in red) 
is the main term 1 — sin(27rr) / (27rr) for the symplectic random matrix 
model for the 1-level density. The lower curve (in blue) includes also 



terms from (13) which are 0(l/log(X)). (In the notation of the Ap- 
pendix, X = 10^^ and AX = 10''.) This version is derived from the 
Explicit Formula. The fit is visibly poor for these values of X and AX. 

As usual, we assumed in ([T]) that supp(^) C {—a, a) C (—1,1), so 
that 

, , / sin(27rr)\ , , , , 

9ir) [-^^) dT = -giO)/2. 

It is really the —g{0)/2 term which appears in the proof via the Explicit 
Formula. Miller [7j derives a version of the 1-level density in which the 



term —g{0)/2 is replaced by a more complicated expression, see (14) 



in the Appendix. This version is shown in green in Figure [2] The fit 



appears to be very good. Since (14) was first derived in [2] from the 



L-functions Ratios Conjecture, the data seems to be good numerical 
evidence for the conjecture. This also seems to indicate is that there is 
a lot of structure in the 0(1/ log(X)) error in ([T]), and the L-functions 
Ratios Conjecture captures that structure. 
The data are available at 

http : //www. math. ucsb . edu/~stopple/quadratic . experiment 



4. Implementation Notes 

Error estimates. With D = 15 digits of accuracy and d near 10^^, 
the crude estimate ^ requires N = 5.4 ■ 10^ terms in the series. We 
can actually do a little better. Using this as a starting estimate, Math- 
ematica^s FindRoot uses a variant of the secant method to determine 
that X = 4.3 ■ 10*^ satisfies for a savings of better than 20%. 

The stronger inequality in the Lemma determines good values for 
the Taylor series truncation parameter B in computing (^(s, vrn^/d). 
Consider the case a = 1, i.e. a negative discriminant. Assuming the 
errors ([T]) in the terms are independent, and using |x('^)| < 1, = n, 
then the standard error in the sum over all in Ij is bounded by [3] 
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We can estimate the sum by an integral 
So the error from the sum over r? in Ij is about 

For d near 10^^, we need T = 65 intervals, and it is easy to compute 




(11) in Mathematica for various B. We see that B = B{j) should 



increase linearly from 84 at j = 31 to 107 at j = 65, in order that the 
total of all errors is only about 10. (For j < 31, the intervals /. 



contain not many more than B{j) squares n^, so the contribution of 
these n, namely 1 < n < 1160, is computed directly.) 

The case a = 0, i.e. positive discriminant is treated similarly. It 
turns out one needs B{j) to increase linearly from 70 at j = 31 to 78 
at j = 65. 

Algorithms. To find fundamental discriminants, I check the congru- 
ence condition, and test for divisibility by the squares of the first 200 
primes. (The 94 examples divisible by the square of a prime larger than 
the 200th prime were easily identified with Mathematica and removed 
from the data by hand.) 

To compute r(s) I use the Lanczos algorithm as in [HIIIZ]- Precom- 
puted values of T{s) allow efficient computation of incomplete Gamma 
functions r(s,a;) for various x by the methods of [9]: series expansion 
for X < 6. and continued fractions for x > 6. These algorithms com- 
pare well with those implemented in Mathematica, giving both absolute 
and relative error no worse than 10. ^^'^ for the relevant range of x and 
|Im(s)| < 1. 

To find zeros of Z{t, x) the computation stepped through values in 
increments of t of size 27r/log(10^^/(27r))/50, i.e. l/50th the mean gap 
between zeros. When a sign change was observed, Ridder's method 
[9] was used to find the root. No effort was made to verify the GRH, 
or that all zeros of Z{t,x) with t < 1. were located. (However, the 
obvious check that Z{0.,x) > 0. was made.) 

Hardware. Computations were done on a 3.0 GHz 8-core Mac Pro. 
Both the integer arithmetic and also the recursion for the derivatives 
G(^)(s,x) were done with GMP 4.2.1 [16j (ported to the Intel Core 2 
Duo 64 bit processor by Jason Worth Martin ^.) For the rest of the 
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floating point computations, the C types long double, long double 
complex sufficed. 



Parallelization. Most of the computation consists of computing the 
values {n^~FjY in ( 10 ). Since this is independent of there is a gain in 
efficiency by computing the quantities Cjk in (10) for 8 discriminants 
at a time. Parallelism is easily implemented using Pthreads. The 
contribution of the intervals Ij is computed in T separate threads for 
8 discriminants at a time. Once all the precomputation is done, the 
zeros of x) for each of the 8 characters x ^i-re computed in 8 separate 
threads. 



Testing. Accuracy of computed zeros was tested three ways: first by 
recomputing well know examples [131 EH moderate sized discrim- 

inants such as —115,147 and —175,990,483. Second, I also implemented 
the method of [T3] directly in Mathematica and compared a few exam- 
ples for discriminants with absolute greater than 10^^, with agreement 
to 15 digits. Third, I compared with the unpublished data from Ru- 
binstein's thesis [in]. This includes 3601 prime discriminants —d with 
10^^ < d < 10^^ + 2 ■ 10^. The data was in agreement with his to the 
10 digits of accuracy he computed. 



5. Appendix: Refined 1 Level Density 

This Appendix closely follows [7| to determine the 1-level density, 
including lower order terms, for the family of quadratic Dirichlet L- 
functions. Instead of considering the set of all fundamental d < X, I 
adapted the proof for 

T{X) = {X < \d\ <X + AX} 

Where Miller treats the case when Xd is an even function, i.e. d > 0, 
I instead considered Xd odd function, —d < 0. Throughout I assumed 
about AX that 

(12) X^/^ log(X) = o(AX) and AX = o(X). 

Refined 1-Level Density (Miller). Let g be an even Schwartz test 
function such that supp(5') C {—a, a), where g denotes the Fourier 
transform of g. Let 
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Then 

(13) 



poo 



logX 



9{r) 



log(7r) + Re 



sin(2 



TXT 



2nr 



tTCT 



r \A logX 



2Re 



2ReA' 



2niT 



+ 



logX 



+ 



dr 
AXi/2 



Of course, to get the OiX^^^hg^ X/AX^^^) error to be o(l/logX) 
we would need to restrict the support of ^ to be C (—1/2, 1/2). 

Figure [3] shows each of the three non-constant terms which are ab- 
sorbed in the 0(1/ log(X)) error in ([T]), all on the same scale of 27rr/ log(X) 
t. The r'/r term (in green) is slowly increasing and very smooth, while 
A' (in blue) is small and wobbly. Observe that when ({1/2 + ij) = 0, 
the contribution at t = ■y/2 of CVC(1 + 2«t) (in red) is positive and 
large. This follows from fTH Theorem 9.6(A)], which says that 

^'^'^ E — + 0(log(t)). 



|i-7l<l 



SO up to a small error, the logarithmic derivative is determined by the 
nearby zeros p. This 'resurgence' of the zeros of ({s) does not play 
much role in the data (t < 1.) presented here. 

Since the fit of the data to even the 'refined' 1-level density is poor, 
we turn instead to the prediction inspired by the L-functions Ratios 
Conjecture. The term — sin(27rr) / (27rr) is replaced by the real part of 



(14) R{t,X) 



tl^(X)logX 



deJ"{x) 



exp 



-2TiiT 



\og{d/-K) 
IokX 



"P / 3 'KIT 

' 4 logX 



C(2)C (i - 



47rzT 
logX 



+ 



c 



Anir 
logX 



(We have simplified the notation from ^ (1-6)]; see also his Lemma 
2.4). Miller shows [3 Lemma 2.1] that on the Riemann Hypothesis, 



giT)RiT,X)dT = -giO)/2 + OiX 
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Figure 3. All three non constant 0(l/log(X)) terms. 
The C'/C term is in red, the T'/T term in green, and the 
A' term in blue. 



and unconditionally with a larger error. Here as usual supp(5') C 
{-a, a). 

In order that the prediction not depend on the specific discriminants 
in J-'(X), we use summation by parts [TJ Remark 2.3] to estimate 



d<X ^ 



logX 



3X fx 



-27rir/log(X) , 

1 n(x^^'^) 

\tt J l-27rzr/log(X) ^ ^ 

and similarly with the sum over d < X + AX. The difference of these, 
divided by ti^(X) = SAX/n^ + 0{X + AX)^/^ jg ^gg^ estimate 



of (14) and denoted -Rest(T, X). Figure |4] shows how _Rest(T, X) (in red) 
compares to — sin(27rr)/(27rr) (blue). 

The graph in green in Figure [2] has — sin(27rr) / {2ttt) replaced by 



RestiX), and also includes the other 0(l/log(X)) terms. 
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